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Abstract 

> 

00 A complete theoretical presentation of the Continuum-mechanical, Anisotropic 

00 Flow model, based on an anisotropic Flow Enhancement factor (CAFFE model) is 

given. The CAFFE model is an application of the theory of mixtures with continuous 



diversity for the case of large polar ice masses in which induced anisotropy occurs. 
(3 The anisotropic response of the polycrystalline ice is described by a generalization 

0^ of Glen's flow law, based on a scalar anisotropic enhancement factor. The enhance- 

^p. ment factor depends on the orientation mass density, which is closely related to the 

^ orientation distribution function and describes the distribution of grain orientations 

^ (fabric). Fabric evolution is governed by the orientation mass balance, which depends 

^ on four distinct effects, interpreted as local rigid body rotation, grain rotation, ro- 

tation recrystallization (polygonization) and grain boundary migration (migration 
recrystallization) , respectively. It is proven that the flow law of the CAFFE model is 
truly anisotropic despite the collinearity between the stress deviator and stretching 
tensors. 
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1 Introduction 



In order to study the mechanical behaviour of large polar ice masses, we use the method 
of continuum mechanics. Generally, ice is treated as an incompressible, isotropic and 
extremely viscous non-Newtonian fluid, and Glen's flow law (Glen 1955, Nye 1952) is 



used as a constitutive equation. However, for thick polar ice masses anisotropic behaviour 
occurs, and therefore. Glen's flow law must be changed. Many efforts have been undertaken 



to deal with this problem (e.g. 


Hutter 


1983, 


Budd and Jacka 


1989 


Jacka and Budd 


1989 


Azuma||1995, Staroszczyk and Morlanc 


112001 


Morland and Staroszczyk||2003 Placidi et al. 


2006 


Gaghardini et al.|2009 


). In this paper we will present a continuum-mechanical model. 



which is ba s ed on earher work by |Faria| ( |2001| |2006a||b| , |Placidi| ( |2004| |2005D , [Placidi 



et al. (2004), Placidi and Hutter ( 2006a[b ). The model is referred to as the Continuum- 



mechanical, Anisotropic Flow model, based on an anisotropic Flow Enhancement factor, 
or "GAFFE model" for short. 

The macroscopic anisotropy of polar ice is due to its microstructure. Ice is a polycrys- 
talline material made of a vast number of crystallites (grains), the mechanical behaviour 
of which is extremely anisotropic ( Jacka and Budd|l989 ). A single ice crystal shows trans- 
versely isotropic behaviour, and its c-axis (optical axis) defines the privileged direction 
(Boehler 1987). Slide along basal planes, orthogonal to the c-axis, is easier than slide 
along other crystallographic planes, and since the study by McGonnel (1891) it has been 
common to refer to this as the deck-of-cards behaviour of ice. However, the transition from 
the mechanics of a single crystal to that of a huge polycrystal entails the complication of 
different deformation mechanisms, and the selection of these mechanisms in different sit- 
uations. A continuum approach is deemed appropriate in order to deal with this problem 
in a manageable fashion. We choose the approach of describing the polycrystalline ice as a 
mixture ( Truesdell||1957a|[b ) , the species of which are the grains characterized by a certain 
orientation (Faria et al. 2006). The orientations of the crystallites, i.e., the unit vectors 
parallel to the c-axes, belong to a continuous space, so that the ice is considered a Mixture 
with Continuous Diversity (MCD) ( |Faria||200lD . 

In the MCD theory, equations are defined at the species level ( "microscopic" ) and at the 
mixture level ("macroscopic"). However, the "microscopic" equations do not govern the 
evolution of single crystallites, the dynamics of which is hidden in the theory. The objective 
of the mixture approach is to predict the polycrystalline behaviour only. In other words, the 
GAFFE model is a macroscopic model, and the microstructure is taken into account only 
phenomenologically, without going down to the actual microscopic level. The same holds 
for classical continuum mechanics: we know that matter is discontinuous, spaces between 
particles (atoms, molecules, etc.) are empty and the molecular structures have strong 
influences on the mechanical behaviour. However, the microstructural characteristics are 
not resolved in detail; instead, constitutive equations for a continuous body are postulated 
that are in accordance with experimental data and general principles like determinism, 
objectivity and the Second Law of Thermodynamics. In the same way, in a MCD, the 
behaviour of a single species is important for the mixture dynamics, but single species 
dynamics does not correspond to any measurable quantity. 

The general set of equations for polycrystalline ice modelled as a MCD was developed 
by [Ma] ( |200l| |2006aD , [Faria et al.| ( |205Bj ), [Ma| ( |2006bD , [Placidi and HutteF] ( |2006b[ ), and 
restrictions of the constitutive equations due to the Second Law of Thermodynamics were 
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given by these authors. Placidi (2004, 2005), Placidi and Hutter (2006a) suggested exphcit 



forms for the constitutive relations. In this study, we present the GAFFE model as an 
improved version of these previous formulations. After defining the notation (Section [2]), 
we derive in Section [3] in a rational way the set of GAFFE equations, while taking care 
that the following requirements are satisfied: 

• All fundamental principles of classical continuum mechanics must be fulfilled. 

• The model must be sufficiently simple to allow numerical implementation in current 
flow models for polar ice masses. 

• The parameters of the model must in principle be measurable by either laboratory 
experiments or field observations. 



In Section 3.1, we define the orientation mass density, and in Section 3.2 we deal with 



the general mass balance that governs its evolution. In Section 3.3, we characterize the 



constitutive quantities introduced in Section |3.2| in order to describe grain rotation, local 
rigid body rotation, grain boundary migration (or migration recrystallization) and polygo- 
nization (or rotation recrystallization) 



In Sections 3.4 and 3.5 



we present an anisotropic 
generalization of Glen's flow law based on a scalar anisotropic enhancement factor. It is 



similar to that by Placidi and Hutter (2006a), but simpler and consistent with the Second 



Law of Thermodynamics. 

Section |4] is devoted to the analysis of the anisotropic properties of the GAFFE flow law. 
We note the coUinearity between the stretching and the stress deviator tensors and show 
that collinearity and isotropy do not share any fundamental concepts, in the sense that 
non-coUinear, isotropic flow laws as well as coUinear, anisotropic flow laws are possible (e.g. 



Liu 2002, Faria 2008). A discussion of the advantages and disadvantages of collinearity in 



an anisotropic flow law terminates Section |4j 

Some examples for the constitutive relations introduced in Section 3^ are provided in 
the appendix (Section [A|). 



2 Notation 

For a general fleld A, the star superscript A* denotes an orientational dependence A* (x, t, n), 
where t is the time, x the position vector and n the orientation (unit vector parallel to the 
c-axis) in the present conflguration. Otherwise the fleld A (x, t) is supposed to be indepen- 
dent of n. It is implicitly assumed that for a given position x all possible orientations n 
are deflned. The gradient (V) and divergence (V-) operators are applied, as usual, to the 
space variable x, while the orientational gradient (V*) and divergence (V*-) operators are 
applied to the orientational variable n. For arbitrary scalars A* and vectors A*, we deflne 



V ■ A* = tr [VA* 



VA* 



tr 



_ dA* 

dA*' 
dx 



V*A* 



V* ■ A* 



dA* 
dn 

= tr 



dA* 








\ dn 


■ nj n 


'dA* 1 


'dA* 


dn \ 


, dn 



n n 



(2) 



3 



An immediate consequence of Eq. ^2 is V*A* - n = 0. We also note that the exphcit form 
of the orientational gradient operator in spherical coordinates is 



ei 



+ 62 



cos 9 cos (f 
cos 9 sin 



dA* sin ip dA* 

09 sin 9 dip 

dA* cos ip dA* 

\ 

d9 sin 9 d(p 



+ 63 



sin 6' 



dA* 
1)9 



(3) 



where {61,62,63} is a fixed orthonormal basis on which we project vectors and tensors, 
and 9 and are the zenith and azimuth angle, respectively. The orientation n can be 
parameterized as follows. 



n 



sin 9 cos (p 
sin 9 sin ip 
cos 9 



sin 9 cos 9? 61 + sin 9 sin 62 + cos 9 63 . 



(4) 



3 CAFFE model 



3.1 Orientation mass density, orientation distribution function 

In the CAFFE model, each point of the continuous body is interpreted as a representative 
volume element of the polycrystal that encloses a large number of crystallites with their 
own orientations. Each orientation is represented by a unit vector n G iS^ (where 
denotes the unit sphere) parallel to the c-axis (Fig. [l^). 




Figure 1: (a) Each point of the unit sphere iS^ represents a particular orientation n. The 
orientation transition rate u* is orthogonal to n. Panel (b) shows a projection on the plane 
spanned by u* and n. The anticlockwise direction of the rotation corresponds to a spin 
velocity s* pointing out of the plane. 

In the MCD framework, distributions of continuous species parameters (like the ori- 
entation) are expressed in terms of associated mass densities. This means that for every 
time t and position x an orientation mass density q* (x, t,n) is defined such that, when 
integrated over 5^, the usual mass density of the polycrystal g (x, t) results, 

0{:>^,t) = J g* {x,t,n) d'^n, (5) 

where d^ra (= sin 9 d9 dcp in spherical coordinates) is the solid angle increment on the 
unit sphere S^. The orientation mass density g*, as stated in Eq. (^, has the following 



4 



physical meaning: the product q* (x, t, n) d^n is the mass fraction of crystalhtes with 
orientations directed towards n within the sohd angle d^n. Therefore, the quantity Q* / Q 
can be interpreted as the analogue of the usual orientation distribution function (ODF) 



(e.g. Rashid 1992, Placidi and Hutter 2004), which is also used in the context of liquid 
crystals (e.g. Blenk et al. 1992, Papenfuss 2000) or in mesoscopic damage mechanics (e.g. 
Massart et al. 2004, Papenfuss and Van|2008 ). However, we remark that in the glaciological 



community, the ODF usually describes the relative number, and not the mass fraction, of 
grains with a certain orientation. 



3.2 Orientation mass balance 

Some kinematic quantities are required in order to describe the evolution of positions and 
orientations. We will use the classical velocity v (x, t) and the orientation transition rate 
u* (x, t, n) as kinematic rates. The velocity v represents the transition of mass from a 
given position to a neighbouring position in three-dimensional space. Analogously, the 
orientation transition rate u* represents the transition of mass from a certain orientation 
to a neighbouring orientation on the unit sphere. Note that the velocity is assumed to be 
independent on the orientation, while the orientation transition rate can depend on the 



orientation [for a longer discussion on this topic see e.g. Faria (2006a) . 

As shown by Faria (2001), the normality condition n ■ n = 1 makes the orientation 
transition rate u* orthogonal to n (u* ■ n = 0) and to the spin velocity s* of the crystallites 
(u* = s* X n); see also the caption of Fig. [T]d. In the MCD theory, the balance of mass is 
formulated as 



dt 



(6) 



The quantity F* is the specific mass production rate. It describes the rate of change of 
mass (per unit mass) of one species (for crystallites having a certain orientation n) into 
another species with a different orientation. This corresponds physically to the effect of 
grain boundary migration (migration recrystallization). 

Integration of Eq. ^ over the unit sphere gives the balance of mass of the poly- 
crystal, 

' " [H = o, (7) 



dt 



provided that Eq. ([s]) and 



/ ^T* d^n = 0, / V* ■ [^*u*l d^n = 



hold. Note that Eq. (|8j)i describes the conservation of mass, and Eq. (|8j)2 is a consequence 
of the Gauss theorem. 



3.3 Constitutive equations for the orientation mass balance 



Placidi and Hutter (2006b) showed that the orientation transition rate u* can be decom- 

+ u: = W ■ n + < , (9) 



posed into two parts. 



u 



u 
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where u*^,^, = W ■ n is the contribution of the local rigid body rotation of the polycrystal, 
W = Skw L is the skew-symmetric part of the velocity gradient L = Vv and u* is a vector 
to be constitutively prescribed. For the latter, we introduce phenomenologically a further 
decomposition, 



(10) 



in order to separate the physical effects of (i) grain rotation ( "deck-of-cards effect", see 
also Fig. [2|, to be modelled by u*j., and (ii) rotation recrystallization (polygonization), to 
be modelled by the orientation flux q* (Godert 2003). 



(a) 



(b) 



/' 









Figure 2: Schematic, two-dimensional representations of representative volume elements. 
Straight arrows inside the volume elements show c-axes of crystallites well oriented for 
deformation. Curved arrows indicate the rotational motion of such crystallites due to 
grain rotation ("deck-of-cards effect"). Arrows outside the volume elements mark the 
state of stretching. In case (a), the state of stretching is pure shear. The c-axes are well 
oriented at angles of 45° with respect to the axes of extension and compression, and the 
c-axes rotate towards the axis of compression and away from the axis of extension. Case 
(b) shows the same configuration rotated anticlockwise by 45° (rotated pure shear) for two 
equivalent volume elements. 



By assuming a linear dependence on the stretching tensor D = SymL, the constitutive 



vector Ugj. takes the form 



u 



.[((D 



n 



n)n — D ■ nl 



111 



where t > is a material parameter [l is called "shape factor" in the theory of rotational 
diffusion; see e.g. Faria (2001)]. Dafalias (2001) noted that, for the special case l = I, 



the unit vector n remains orthogonal to the associated material area element, and thus 
the rotation is an affine rotation. However, the possibility of non-affine rotations cannot a 
priori be ruled out. An advantage of the present macroscopic approach is the possibility 
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to parameterize t without any conflicts with "microscopic" assumptions. [As a contrasting 



example, the model by Staroszczyk and Morland (2001) is also a macroscopic model, but 



it is restricted to affine rotations.] In fact, Placidi (2004) showed that the fabrics in the 



upper 2000 m of the GRIP ice core in central Greenland can be best explained by the 
value L ^ 0.4. A study on the EPICA ice core in Dronning Maud Land, East Antarctica, 
provided a best fit between modelled and measured fabrics for l = 0.6 (Seddik et al. 2008). 



Larson 1988, Blenk et al. 


1992, 


Larson 


1999 


Papenfuss 


2000, 


Daf alias 


2001 



to ice (e.g. 

and references therein). They were derived by Placidi (2004) within the MCD framework. 
We remark that, even though the unit vector n specifying the orientation of the crystals is 



unique, Eq. (11) is not the most general case. In fact, Dafalias (2001) discussed the case of 



non-affine rotations. More generally, Faria (2001) and Placidi and Hutter (2006b) derived 



the thermodynamically consistent class of these constitutive equations. 



Following the argumentation by Godert (2003), the orientation flux (which is supposed 



to describe rotation recrystallization) is modelled as a diffusive process. 



;i2) 



where the parameter A > is the orientation diffusivity, and T-L* is an orientation-dependent 
"hardness" function. However, recent results by Durand et al. (2008) suggest that rotation 
recrystallization is an isotropic process not affected by the orientation. In this case, the 
choice 

n* = 1 (13) 



is indicated, which renders Eq. (12) equivalent to Pick's laws of diffusion on the unit sphere. 
We remark that in the MCD theory the hardness function is called the chemical potential 
for the given species. It is a constitutive quantity that depends at least on the OMD g*. 
Consequently, it is possible to show by applying the chain rule that often an equivalent of 
Pick's law results even when "H* is not a constant (e.g. Faria 2001). 

As for the specific mass production rate P*, in the studies by Placidi (2004) and Placidi 



(2005) it was shown that, in order to model the effect of grain boundary migration, a 



reasonable constitutive equation for P* is 

F* = f [D* - {D*)] . (14) 
The dimensionless quantity D* is called the stretching deformability of crystallites, 

(D ■ n)2 - ((D ■ n) 



D* 



n 



tr (D2) 



(15) 



its physical meaning is the square of the resolved shear strain rate (or stretching) on 
the basal plane, normalized by the orientation-independent scalar invariant tr(D2). The 
additional factor 5 is merely a convention, for which the reason will become clear below 



(Section 3.4). Further, (■) is the averaging operator 



Q 



(16) 
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and r is a constitutive parameter (see Fig. [3j). The conservation of mass expressed by 
Eq. ([8])i is compatible with Eqs. (14)-( 16), and they are also co mpatible with the ratio- 



nal constitutive theory developed by 



Placidi and Hutter 



(2006b). Provided that T > 0, 



Eqs. ([I4j)-(|16]) have the effect that T 
itjD 



is greater than zero when the stretching deformabil- 
is high, and F* is less than zero when the stretching deformability D* is low. This 

e.g. |Kamb||1972D. 



means that crystallites well-oriented for deformation will be enlarged 
Since ice crystallites deform essentially by basal shearing, the resolved shear rate (which 
is proportional to the orientation dependence of D*) is related to the rate of accumulation 
of deformation energy in the material, which drives dynamic recrystallization. 




Figure 3: Specific mass production rate F* according to Eq. (14) 



The first contribution in Eq. ([9]), u*^^^., is thermodynamically reversible, because there 
is no energy dissipation associated with local rigid body rotations. The second contribu- 
tion, u*, has been split up in Eq. (10) into u*j, and q*. The grain-rotation part, u*j,, is 



thermodynamically irreversible because it is linearly dependent on the stretching tensor 

D, see Eq.JTTL ^ 

fluid (e.g., |Hutter||1983l |Muller||1985l iFariapOOll). 



which by definition must vanish for any reversible process in a viscous 

The rotation-recrystallization part, q*. 



is thermodynamically irreversible because of the diffusive nature of the process as specified 
in Eq. (12). The specific mass production rate F* is also thermodynamically irreversible. 



because grain growth and recrystallization are thermally activated, irreversible processes 



[see also the discussion by Faria et al. (2006, Sect. 3c)]. Note that, besides the physical 
interpretation, the thermodynamic reversibility or irreversibility of these terms can also be 



investigated by exploiting the Second Law of Thermodynamics (Faria 2006a, Placidi and 



Hutter 2006b). 



A problem is that it is not possible at the moment to constrain the values of the two 
parameters F and A in a reasonable fashion. Determining these parameters by experiments 
is very difficult, because the relevant time scales are far too large and the strain rates far 
too low to be reproduced in the laboratory. Deformation experiments, even if conducted 
over a period of years, inevitably end up by activating non-natural deformation and re- 
crystallization mechanisms. The only promising way out of this is to measure the fabrics 
(c-axis distributions) and the changes in grain stereology (sizes and shapes) in natural 
polar ice and fit F and A to these observations. The situation is complicated further by the 
fact that a functional dependence of these parameters on temperature and/or dislocation 
density should be considered (cf. [Faria 2006a|[b ). This requires further attention. 

Some examples for the orientation transition rate (grain rotation) and the orientation 
production rate (recrystallization) under different deformation regimes are given in the 
appendix (Section [A]). 
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3.4 Anisotropic flow law 



The anisotropic flow law presented by Placidi and Hutter (2006a) has been modified in 



order to make it simpler and compatible with the Second Law of Thermodynamics, 



D = AE{S) 




(n-l)/2 



S, 



(17) 



where A and n are the same rate factor and stress exponent, respectively, as in the isotropic 
Glen flow law, S is the stress deviator deflned by 

(t)' ,18) 

(i.e., the deviatoric part of the Cauchy stress tensor t) and I is the identity tensor. Further, 
S E [0, 5/2] is the positive scalar 



S= (S*) = I —S*d?n, 



(19) 



and S* is the analogue of D* 



S* = b 



(S-n)^-((S-n)-n) 
tr(S2) 



(20) 



which can also be written in terms of the Cauchy stress tensor as 



5"*-''!:">'>0. 



tr 



(21) 



We call S* the stress deformability of crystallites and S the stress deformability of the 
polycrystal. In the literature, the scalar (t ■ n)^ — ((t ■ n) ■ n)^ = (S ■ n)^ — ((S ■ n) ■ n)^ 
has been identifled with the square of the resolved stress on the basal plane, so that the 



stress deformability of crystallites, Eq. (20), can also be called the normalized square of 



the resolved stress on the basal plane. As for the stress exponent, it is often chosen as 
3 (e.g. Paterson 1994), but we will keep it general in the following. 



n 



For a thermodynamicist it may appear strange to formulate a constitutive equation in 
terms of the stretching tensor and not in terms of the stress deviator. However, there is no 
inconsistency in this formulation. From the theoretical point of view the stress is indeed 
the constitutive property and the strain rate (stretching) is the variable, but whether we 
choose this or the inverse relation is just a matter of taste or custom. In the glaciological 



community the inverse form is most commonly used; in the book by Hutter (1983) the 
historical reason for this is given. 

Our new formulation of the flow law is not only compatible with the Second Law of 



Thermodynamics (Placidi and Hutter 2006b), but Eq. (17) is much more flexible than 
the previous Placidi-Hutter formulation. The mechanical non- linearity (tr(S2)/2)^" 
and the anisotropic part E{S) are now nicely separated, so that the choice of the stress 
exponent is not limited to n = 3 any more, and the new formulation is not even restricted 
to a power law. 
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In Section 3.5 we will show that, due to Eq. (17), the two quantities D* and S* defined 



in Eqs. (15) and (20) are identical 



S* = D\ (22) 

so that we will simply call them the species (or "crystallite") deformability. In the same 
way, the positive scalars S and D will be called the polycrystal deformabihty, 

S = {S*) = {D*) = D . (23) 

Both the crystal and polycrystal deformabilities can only assume values in the range from 
zero to 5/2, 

5* = G [0, |] , S = De [0, |] . (24) 

The proof (which is laborious and shall not be detailed here) involves to insert Eq. (|4]) in 
Eqs. (15) and (16), and study the maxima and minima of the deformabilities D* and D 



for general stretching tensors D as functions of the zenith angle 9 and the azimuth angle 
0. Taking into account that for a randomly distributed OMD (isotropic fabric) 



Q 



S = D 



(25) 



holds, the function E{S) is demanded to be monotone, of class C^[0, 5/2] and has the 
fixed points 



E(0) = E^in, ^(1) = 1, ^(§) = ^^ax, (26) 

where -Emin < 1 and -Emax > 1 are the minimum and the maximum enhancement factors. 
This means that if the polycrystal deformability S is highest (5* = 5/2), the fiow law (17) 
gives the maximum stretching, if the polycrystal deformability S is lowest {S = 0), the 



fiow law (17) gives the minimum stretching, and if the polycrystal deformability is the 



same as for the isotropic case {S = 1), the fiow law (17) reproduces the classical Glen fiow 
law. 

As for the detailed functional form of the anisotropic enhancement factor E{S), some 
experimental data suggest that the enhancement factor depends on the "averaged Schmid 
factor" to the fourth power ( Azuma||1995 Miyamoto|1999 ) . Since the polycrystal deforma- 
bility S of the GAFFE model is related to the square of the averaged Schmidt factor, it 
is reasonable to assume a dependency of E on S"^. However, this does not allow to fulfill 
Eq. (26 ) for arbitrary choices of the parameters £^mm and -Emax- Hence the function E{S) is 



chosen to depend on S in the interval [1, 5/2] [in which the experiments by Azuma (1995) 



and Miyamoto (1999) have been carried out] only, and for the interval [0, 1] a dependency 
on S'^ is introduced. The exponent r is adjusted such that the function is continuously 
differentiable at S" = 1. This yields 

' E — r 

21 

AS^ (E^ax - 1) + 25 - 4E^ax 

21 



1 EjYiin )S^ + E 

n 



r 



E{S) 



E^ 



[0,1] , 



(27) 



S G 



1 5 



Russell-Head and Budd 


1979 


Pimienta et al. 


1987, 


Budd and 



ndicate that the parameter -Emax (maximum softening) is ~ 10. The parameter 
E'min (maximum hardening) can be realistically chosen between and ~ 0.1, a non-zero 
value serving mainly the purpose of avoiding numerical problems. The function (27) is 
shown in Fig. |4} 
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LU 



E =10 




max 




E . =0.1 




mm 





0.5 



1 1.5 
Deformability 



2.5 



Figure 4: Anisotropic enhancement factor E{S) as a function of the deformabihty S ac- 
cording to Eq. (27), for -Emax = 10 and -Emin = 0.1. 



3.5 Inversion of the anisotropic flow law 



The anisotropic flow law (17) can be inverted analytically as long as the power-law form 
is retained. From Eq. ((TtI) we find 



A' E\S) 



and thus 



tr(D2) Ar(S2)\""' tr(S2) 

2 \ \ V" 



{2t 



(29) 



Inserting this result in Eq. (17) and solving for S yields 



S = A-i/"E-i/"(5) 



tr(D^ 



^(l-l/n)/2 



D. 



(30) 



In order to complete the inversion, we must prove Eqs. (22) and (23). By writing 
Eq. (17) in short as D = 8/(2?]) (where rj is the shear viscosity of the flow law), we obtain 
for the crystalline deformability S* 



S* 



(S ■ n)^ - ((S ■ n) ■ n)^ 

tr(S2) 

{2r]'D ■ n)2 - {{2r]'D ■ n) ■ n)^ 
tr[(2r/D)2] 

(D ■ n)^ - ((D ■ n) ■ n)^ 
tr(D2) 



(31) 



so that Eq. (22) is proven. By applying the averaging operator (16) to this result, we find 
immediately 

S = D, (32) 
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which is the assertion of Eq. (23). Hence, we can replace S* by D in Eq. (30), which yields 
the inverted anisotropic flow law 



S = A-i/"E-^/"(^) 



tr(D2 



-{l-l/n)/2 



D. 



(33) 



4 On the anisotropy of the CAFFE flow law 



In Sections 3.4 and |3.5| we have given the generalization of Glen's flow law considered in 
the CAFFE model. The anisotropy of the CAFFE flow law has recently been analyzed 



mathematically by Faria (2008) through the derivation of the symmetry group of the 
CAFFE model. In this section, we prove in a more direct way that the CAFFE flow is 
anisotropic despite the collinearity between the tensors S and D, and give examples in 
order to justify this choice. 

4.1 CAFFE anisotropy in the context of material theory 

In the context of constitutive theory, the definition of isotropy states that any rotation 
of the body in question does not alter its material response. Mathematically speaking, 
this means invariance of the material functions (or functionals) to arbitrary orthogonal 

86), so that the 



Liu 2002 



P- 



transformations P of an undistorted configuration k (e.g. 
symmetry group Q of the material contains the entire group of orthogonal transformations 
0(3) [Q ^ 0(3)]. Anisotropy is the logical opposite: for at least one orthogonal transfor- 
mation the invariance does not hold, so that the symmetry group does not contain the 
entire group of orthogonal transformations [Q ^ 0(3)]. 



By construction, the anisotropy of the CAFFE flow law (17) must be contained in the 



enhancement factor E{S) via the polycrystal deformability S. So let us assume that, at 
the initial time t = 0, the initial configuration Kt=o is given by an unloaded ice specimen 
with the OMD f)*(x, n, 0). At t = 0"*", it is subjected to the stress S, and, according to 



Eqs. (19) and (20), the resulting polycrystal deformability is 



S 



f)tr(S2) Js 



g [n) 



:s-n) 



((S ■ n) 



n 



d'n. 



(34) 



where, for simplicity of notation and in the rest of this session, we will omit the dependence 
of OMD on position and time. Now let us consider a second initial configuration Kt=o 
rotated by a proper orthogonal tensor P with respect to ^^=0- The rotated orientations 
are given by 

n = P ■ n (35) 
(Fig. [sj. The OMD follows the rotation, so that 



(nj 



>T-fi) 



(36) 



At t = 0"^, the rotated configuration is subjected to the stress S, which is supposed to be 
the same as before, 

S = S (37) 
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(Fig. [5]). Tlie polycrystal deformability with respect to the rotated configuration is then 

5 



S 



te4 



pell, pTt 



^tr(S2 



^tr(S2 



^*(n) (S ■ fi)^ - ((S ■ fi) ■ fi)2 



{31 



Let us change the name of the integration variable in the last integral of Eq. (38) from n 
to n, 

5 



S 



^tr(S2 



g*{P' -n) (S-n)'-((S-n)-n) 



d n . 



(39) 



This is the same as the polycrystal deformability with respect to Kt=o [Eq. (|34j)] for arbitrary 
transformations P G 0(3) if and only if ^*(n) = const = g/^An). In this case, the fiow 



law (17) is isotropic. For the general case of a non-constant OMD, the deformabilities 



(34) and (39) are not equal for arbitrary transformations P, so that the fiow law (17) is 



anisotropic, QED. From a mathematical point of view, Eqs. (34) and (39) make clear that 



the symmetry group Q of the material defined by the anisotropic GAFFE fiow law includes 
the invariance group of orthogonal transformations that keep the orientation mass density 



g* unchanged (e.g. Faria 2008). 




Stress S 



Response S 



4)^ 



Stress S = S 



("same experiment") 
n = P • n 



Response S 



Figure 5: Anisotropy of the GAFFE fiow law: If the same stress (S = S) is applied to two 
rotated initial configurations {Kt=o, iit=o), the responses S and S are different in general. 



4.2 Anisotropic behavior for simple shear 

Let us illustrate the anisotropic behavior of the GAFFE flow law by a simple example. For 
the vertical single maximum fabric 

^*(x,n,t) = f?(x,t)5(n-e3) (40) 

and the simple shear deformation 

2 

D = I I , (41) 

1 
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we find 



6{n - 63) 




S 



(42) 



Thus, the x-z component of the fiow law (17) yields 

7 = 2AE(|)r" = 2AE„ 



(43) 



where the explicit definition of the function E{S) [Eq. (27)] has been used. 

Now we rotate the sample around the direction 62 by 45° and keep the experimental 
apparatus fixed [apply the same stress in the sense of Eq. (37)]. In this case, the OMD 
changes as follows, 

ei3 = 

V A J 



g* = g S{n - eis). 



(44) 



while the state of stress (42)2 is unchanged. The species deformability is equal to zero for 
the orientation 613, 

S*{eis) = S = 0. (45) 



It follows that the stretching tensor evaluated with the fiow law (17) yields the shear rate 

7 = 2AE{0) r" = 2AE,^i^T'' . (46) 
Since £"111111 ^ -E'max, the shear rate of Eq. (46) is much smaller than that of Eq. (|43|). In 



other words, the material response of the ice specimen has changed considerably due to the 
45° rotation of its initial configuration. This fulfills clearly the criterion for an anisotropic 
material. 



4.3 Isotropic, anisotropic, collinear and non-coUinear flow laws 

On the one hand, the classical Glen fiow law 



D = A 



(n-l)/2 



s, 



(47) 



which results from the GAFFE fiow law (17) by setting E{S) = 1, is isotropic and collinear 
with respect to the tensors D and S. On the other hand, many anisotropic fiow laws pub- 



lished so far relate D and S by tensor quantities {. 


Lliboutry||1993 


Azuma 


1994 


Mangeney 


et al. 11996, 


Svendsen and Hutter||1996 


Godert and Hutter||1998 


rhorsteinsson 


2001, 


Mor- 


land and Staroszczyk 


2003 


Gillet-Ghaulet et al. 


2005 


), thus giving up the collinearity 



between D and S. 

This often leads to the misconception (at least in the glaciological community) that 
isotropic fiow laws must be collinear and anisotropic fiow laws must be non-coUinear with 
respect to D and S. However, this is not the case. As we have seen above, the GAFFE 
fiow law is anisotropic, but collinear. Gonversely, a non-collinear fiow law is not necessarily 
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anisotropic. An example is the general Reiner-Rivlin flow law for isotropic viscous fluids 
(e.g. Liu 2002, p. 109). For the incompressible case, it reads 



S = aiD + a2 [B^ ij , (48) 

where ai and 0^2 are material parameters. Provided that 02 7^ 0, this flow law is evidently 
non-coUinear. So we highlight that isotropy/anisotropy and collinearity/non-coUinearity 
are two entirely different concepts, and all four possible combinations can be realized. 

The disadvantage of using a collinear anisotropic flow law is that, for a given fabric 
and a given state of stretching, we select a single viscosity of the polycrystal that is the 
same for every components of the stress deviator. However, in reality, for complex states 
of stretching (superposition of compression and shear etc.) different directions will show 
different degrees of softening or hardening. This shortcoming is a tribute to the simple 
formulation with a scalar enhancement factor, which allows to set up the flow law with 
only two well-known parameters (-Emin, -f'max)- 



5 Conclusion 



We have presented a constitutive model for the dynamics of large polar ice masses. This 
GAFFE model consists of an anisotropic generalization of Glen's flow law based on a scalar 
enhancement factor, and a fabric evolution equation based on an orientation mass balance. 
The latter arises from the framework of Mixtures with Gontinuous Diversity and uses the 
orientation mass density as the variable which describes the anisotropic fabric. Three 
constitutive quantities have been introduced, namely the orientation transition rate due to 
grain rotation, the orientation flux and the specific mass production rate. They have been 
linked to the physical processes of grain rotation, rotation recrystallization (polygonization) 
and grain boundary migration (migration recrystallization), respectively. The anisotropy 
of the GAFFE flow law has been proven, and in that context it has been emphasized that 
isotropy/anisotropy and collinearity/non-collinearity (between the stress and stretching 
tensors) must be clearly separated. Some applications of the GAFFE model to simple 
deformation states have been discussed (in the appendix), for which analytical solutions 
could be obtained, and which could easily be checked for their physical plausibility and 
consistency with observations. 

Due to its relative simplicity, the GAFFE model is suitable for implementation in ice- 



flow models. This has already been done by Seddik et al. (2008) for a one-dimensional 



model of the site of the EPIGA ice core at Kohnen Station in Dronning Maud Land, East 



Antarctica (EPIGA Gommunity Members 2006), and by Seddik (2008) and Seddik et al. 



(2009) for the three-dimensional, full-Stokes model Elmer/Ice in order to simulate the ice 



flow in the vicinity within 100 km around the Dome Fuji drill site (Motoyama 2007) in 
central East Antarctica. 
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A Examples for the evolution of the orientation mass 
density 

A.l Evolution due to grain rotation 

The deck-of-cards deformation mechanism (grain rotation) implies that the c-axis of a crystallite 
in the polycrystalline aggregate rotates towards the axes of compression and away from that of 
extension. This is illustrated graphically in Fig. [2] for the case of rotated pure shear (simple 
shear) and described mathematically by Eq. ( |11[ ), provided that the constitutive parameter l is 
positive. Equation (11) fulfills the principle of material frame indifference. Thus, if the rules 



for compression and for extension are satisfied, then the rules for simple shear are a direct 
consequence. Here we give some simple examples in which this can explicitly be seen. 

We use a Cartesian frame of reference for which the orientation n of the c-axis is parameterized 
by Eq. Q. For uniaxial vertical compression (transversely isotropic horizontal extension), the 
stretching tensor is 







D 



2 

- 









1 1 

- e ei ei - e e2 e2 



(49) 



where e > holds. From Eqs. Q, (11) and (49) we derive the explicit form of the orientation 
transition rate due to grain rotation. 



u 



gr 



COS (p COS 6 

i£sm.29 I sin 99 cos 



sm ( 



L £ sin 29 (cos ip cos ^ ei -|- sin (p cos 9e2 — sin 9 63) . 



(50) 



The direction of u*j, is coherent with the rules of Fig. |2p, (see also Fig. |lp): The third component 
Ugj, • e3 is positive when 9 < 7r/2 and negative when 9 > 7r/2. If the crystallites are in the plane 



spanned by ei and {(f = 0), Eq. (50) simplifies to 

. 3 



/ — cos 9 



V9 = 



u 



gr 



i £ sin 29 



\ 





sin( 



(51) 



and if they are in the plane spanned by 62 and 63 {ip = tt/2), we find 

3 ^ 

p = — =^ u*^ = - i£ sin 29 



TT 

2 



V 




(52) 
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It is worth noting that the norm of u*j. which results from Eqs. (50), (51) or (52) shows the 
exphcit dependence on the Schmidt factor sin 20 that Azuma and Higashi (1985^) recognized 
experimentally. The presence of the Schmidt factor guarantees that crystallites with vertical 
and horizontal orientations do not rotate, while those at orientations 45° off the vertical show 
maximum rotation. 

If the uniaxial compression is along the first (ei) instead of the third (e^) axis of the frame 
of reference, the orientation transition rate due to grain rotation takes the form 

/ 



D 



-e 



V 



£ 

2 







E 

2 



gr 



sin 9 cos (f + B {9, ip) sin 9 cos 
te \ — 2 sin 9 sin (p + B {9, (p) sin 9 sin (p 
-\cos9 + B {9,(p) cos 6* 



where 



B{9,^) 



sin 9 cos p sin 9 cos ip + sin 9 sin p sin 9 sin p + cos 6 cos ( 



(53) 



(54) 



The difference to Eq. (50) arises only because the spherical coordinate system, which underlies 



the representation of the unit vector n in Eq. (Q , is more convenient for the vertical compression 
(49) than for the horizontal compression (|53|)i. 

For a planar elongation (or pure shear) state of deformation, 



D 









e ei ei - e es es 



(55) 



the orientation transition rate due to grain rotation results from Eqs. (|4]), (11) and (55) as 

(— 1 + sin 9 cos p> sin 9 cos p> — cos 9 cos 9) sin 9 cos ip 



u„ 



\ 



- 1 + sin 9 cos p sin 9 cos ip — cos 9 cc 


sin 9 sin 0(1 + cos p? cos ip) cos 9 



(56) 



which, in the plane spanned by ei and 63, is 



gr 



i £ sin 29 




(57) 



This is larger by the factor 4/3 compared to the orientation transition rate of Eq. (51)2- The 



reason for this difference is that the component D22 does not contribute to grain rotation in the 

1 by ei and 63, which makes the orientation trans 
(where D22 = 0) faster than in the case of Eq. (51) (where D22 



plane spanned by ei and 63, which makes the orientation transition rate in the case of Eq. (57) 

_ ^/2). 
For the simple shear situation of Eq. (41), which is illustrated in Fig. 12" 



we find for the 



orientation transition rate due to grain rotation 



u„ 



7 



/ COS0 (2 sin^ ^cos^ 99 — 1) 
2 sin 9 sin 29 sin 2ip 



(58) 



y sin 9 cos p cos 29 
which, in the plane spanned by ei and 63, gives 

/ - 







L ^ COS 29 



\ 



- cos I 


sin 9 



(59) 



The direction of u*j. is once more consistent with the rules of Fig. [2|3 (see also Fig. [Tj?). For 
instance, for crystallites oriented upward within 9 < 7r/4 and a positive shear rate 7 > 0, the 
component of Ugj, along ei is negative. 
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A. 2 Evolution due to recrystallization 



We now give some examples for the recrystallization term of Eq. (14) for standard deformation 
situations, and also provide a model of the experiments by Budd and Jacka ( 1989 ) for uniaxial 
compression with isotropic horizontal extension. 

For the case of pure shear rate as defined in Eq. (55), crystallites oriented vertically have a 
vanishing species deformability. 



n 



63 



D* 



(D • n)2 - ((D • n) • n) 



3p2 

1 ^ 



0. 



tr(D2) 2 

while crystallites inclined by 45° off the vertical towards ei have the maximum deformability. 



(60) 



n 



ei3 




D* 



(D • n)2 - ((D • n) • n)^ 







(61) 



tr(D2) 2e2 

For general anisotropic fabrics, the averaged deformability {D*) is between these extremes. It 
follows from Eq. (14) that the favourably oriented crystals with n = eis will grow (F* > 0) and 
the unfavourably oriented ones with n = es will shrink (F* < 0). This is the physically expected 
behaviour. 

The experiments by |Budd and Jacka (1989) were carried out under the deformation regime 
of uniaxial compression with isotropic horizontal extension, as specified by Eq. (49). By using 
the parameterization (j4| for general orientations n, we compute the species deformability (15) as 

15 



D* 



■ sir? 6 cos^ ( 



If at the initial time t = 0, the OMD is random {q* = q/^'k, isotropic fabric), then 



D = {D*) 



— Dan 







27r 2 

47r 



dip 



15 



sin2(9cos2 6' sin6'd(9 = 1 



(62) 



(63) 



where d^n = sin 9 dO dip and usual integration rules have been used. The specific mass production 
n Eqs. 

T*{e) = F 



rate which results from Eqs. (14) and (63) is 

15 



sin^ 6 cos^ I 



15 



8 



f sin^ 26 -V . 



(64) 



2 arcsm ^, 



grow \r*{e) > 0], while the 



Consequently, crystallites with 6 & [ \ arcsin y^, f 
others shrink [F*(0) < 0], and an anisotropic fabric evolves. If we do not consider grain rotation 
and rotation recrystallization, then, asymptotically for t — )• oo, we will obtain 



2tt sin ^1 



'0 



5(0 -^o), with 00 = 45° 



(65) 



where 5 is the Dirac delta function. Equation (65) is the mathematical representation of a girdle 
fabric (see e.g. Placidi and Hutter 2006a) in which all the crystallites are inclined by 45° with 
respect to the vertical. If grain rotation is superimposed, these crystallites will experience an 
additional rotation towards the compression axis 63 [in accordance with Eq. (59) or Fig. [2^], so 
that the small girdle fabric observed by Jacka and Budd ( 1989 ) is deduced. 
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Another interesting example is the rotated pure shear (simple shear) regime of Eq. (41) for 
general orientations n represented by Eq. ([4]). The crystal deformability is computed for this case 



as 



D* 



(66) 



cos (/7sin 6 + cos 0(1 — 4 cos 99 sin 
If at the initial time t = 0, the OMD is random {g* = q/Att, isotropic fabric), then, analogue to 

D = 1 (67) 

5 , 



Eqs. (63) and (64), we find 
and 



cos ip sin 9 + cos 6* (1 — 4 cos ips'm i 



(68) 



For times t > 0, an anisotropic fabric evolves, because crystallites oriented near n = ei (9 ^ ir/2 
and 93 ~ 0) or n = 63 (0 ~ 0) grow and the others shrink. Hence, without local rigid body 
rotation, grain rotation and rotation recrystallization, asymptotically for f — t- 00 we will obtain 
the two-maxima fabric 

g* = ^6{n-ei) + ^6{n-e3) . (69) 



If grain rotation and rigid body rotation are superimposed, the fabric reported by Kamb ( 1972 ) 
results. 
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